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1. Introduction 

Hydromagnetic processes play an impor- 
tant role in many astrophysical systems 
(e.g. stars, galaxies, accretion discs). This 
is because the medium is hot enough to be 
partially or fully ionized. Because of the 
huge scales involved the medium is usually 
turbulent, provided there is an instability 
(shear, convection) facilitating the cascad- 
ing of energy down to small scales. 

In turbulence research it has been a long 
standing tradition to solve the hydrody- 
namic equations using spectral schemes 
which have the lowest possible discretiza- 
tion error. Spectral schemes are particularly 
useful for incompressible problems where 
one needs to solve a Poisson-type equation 
for the pressure. However, spectral schemes 
are no longer optimal in many astrophysi- 
cal circumstances where flows are generally 



compressible. Lower order spatial deriva- 
tive schemes are generally unacceptable 
in view of their low overall accuracy, even 
when schemes are used where mass, mo- 
mentum, and energy are conserved to ma- 
chine accuracy. On massively parallel ma- 
chines, on the other hand, spectral schemes 
are difficult to make run efficiently. High or- 
der finite difference schemes are therefore a 
useful compromise. Such schemes can yield 
almost spectral-like accuracy. 

Our code uses centered finite differences 
which make the adaptation to other prob- 
lems simple. Since the code is not written 
in conservative form, conservation of mass, 
energy and momentum can be used to mon- 
itor the quality of the solution. A third or- 
der Runge-Kutta scheme with 2N storage 
is used for calculating the time advance. 



2. Advantages of high-order schemes 

Spectral methods are commonly used in 
almost all studies of ordinary (usually in- 
compressible) turbulence. The use of this 
method is justified mainly by the high nu- 
merical accuracy of spectral schemes. Al- 
ternatively, one may use high-order finite 
differences that are faster to compute and 
that can possess almost spectral accuracy. 
In astrophysics, high-order compact finite 
differences have been used to model stel- 
lar convection and shear fiows in ac- 
cretion discs 0. In contrast to explicit fi- 
nite differences, compact finite differences 
have a smaller coefficient in the leading 
error, even if both schemes are of the same 
order. However, compact schemes are still 
nonlocal in the sense that each point affects 
every other point, which enhances commu- 



nication. This is the main reason why we 
adopt explicit centered finite differences. 

In this section we demonstrate, using sim- 
ple test problems, some of the advantages 
of high-order schemes. The explicit formu- 
lae for first and second derivatives are 



+45/,+i-9/,+2 + /.+3)/(605x). 



= (2/,_3 - 27/,_2 + 270/,„i - 490/, 
+270/,+i - 27/^+2 + 2/,+3)/(1805x2). (2) 

Full details of these schemes, including for- 
mulae for the boundaries, can be found in 
Ref. 0. This scheme was also used in re- 
cent applications to the problem of resis- 
tively limited growth in models of stellar dy- 
namos; see Ref. [0 . 

It is commonly believed that high-order 
schemes lead to Gibbs phenomena and that 
more viscosity is needed to damp them out. 
In fact, the opposite is true as is demon- 
strated in Fig. |1| for advection tests and in 
Fig. for the stationary Burgers shock. 

For the time stepping high-order schemes 
are necessary in order to reduce the am- 
plitude error of the scheme and to allow 
longer time steps. Usually such schemes re- 
quire large amounts of memory. However, 
there are the memory-effective 2A^-schemes 
that require only two sets of variables to 
be held in memory. Such schemes work 
for arbitrarily high order, although not all 
Runge-Kutta schemes can be written as 
2A^-schemes These schemes work it- 
eratively according to the formula 

Wi = aiWi-i + 6tF{ti_i,Ui-i), (3) 

Ui = Ui_i + (3,iWi. (4) 
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Fig. 1. Advection tests with schemes of different spatial or- 
der. Resulting profile after advecting a step-like function 5 
times through the periodic mesh. The dots on the solid line 
give the location of the function values at the computed 
meshpoints and the dotted line gives the original profile. 
For the panels on the right hand side the diffusion coeffi- 
cient is too small and the profile shows noticeable wiggles. 
5x = 1/60. 

For a three-step scheme we have i = 1, 3. 
In order to advance the variable u from u^"^ 
at time t^") to at time t("+^) = t^") + 

Shwe set in Eq. (H) 

and = Us, (5) 



Mo 



with Ui and U2 being intermediate steps. In 
order to be able to calculate the first step, 
i = 1, for which no Wi-i = wq exists, we 
have to require cti = 0. Thus, we are left 
with 5 unknowns, 02, «3, Pi, P2, and P3. 
Three conditions follow from the fact that 
the scheme be third order, so we have to 
have two more conditions. One possibility 
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marginally better in some simple test prob- 
lems where an analytic solution was known. 

Table 1 

Possible coefficients for different 2Af-RK3 schemes. 



label 



0-2 



Pi /32 /33 



symmetric -2/3 -1 1/3 1 1/2 

predictor/corrector -1/4 -4/3 1/2 2/3 1/2 
inhomogeneous -17/32 -32/27 1/4 8/9 3/4 
Williamson (1980) -5/9 -153/128 1/3 15/16 8/15 
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Fig. 2. Burgers shock with schemes of different spatial 
order and different values of the viscosity. The solid lines 
give the analytic solution. A second order scheme (top 
row) requires a viscosity of at least O.fiuSx, where u is 
the amplitude of the shock and 6x the mesh spacing. For 
a sixth-order scheme, a viscosity of O.fiuSx yields good 
results, and for a tenth-order scheme a viscosity of O.bu&x 
can be used. 

is the choose the fractional times at which 
the right hand side is evaluated, for example 
(0, 1/3, 2/3) or even (0, 1/2, 1). In the latter 
case the right hand side is evaluated twice 
at the same time. It is therefore some sort 
of predictor-corrector scheme. Yet another 
possibility is to require that inhomogeneous 
equations of the form -it = with n = 1 and 
2 are solved exactly. The corresponding co- 
efficients are listed in Table |l| and compared 
with those given by Williamson |jl| . In prac- 
tice all of them are about equally good when 
it comes to real applications, although we 
found the first one in Table |I] ('symmetric') 



3. Implementing magnetic fields 

Implementing magnetic fields is relatively 
straightforward. On the one hand, the mag- 
netic field causes a Lorentz force, J x B, 
where B is the fiux density, J = V x B//io 
is the current density, and /xq is the vacuum 
permeability. On the other hand, B itself 
evolves according to the Faraday equation. 



(9B 
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-V X E 



(6) 



where the electric field E can be expressed 
in terms of J using Ohm's law in the labo- 
ratory frame, E = — u x B + J/a, where 
a = (?7/io)~^ is the electric conductivity and 
rj is the magnetic diffusivity. 

In addition we have to satisfy the condi- 
tion V ■ B = 0. This is most easily done by 
solving not for B, but instead for the mag- 
netic vector potential A, where B = V x A. 
The evolution of A is governed by the un- 
curled form of Eq. (|^) , 

dA 
'dt 

where is the electrostatic potential, which 
takes the role of an integration constant 



-E- Vc 



(7) 
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which does not affect the evolution of B. 

The choice = is most advantageous 
on numerical grounds. (By contrast, the 
Coulomb gauge V ■ A = 0, which is very 
popular in analytic considerations, would 
obviously be of no advantages, since one still 
has the problem of solving a the solenoidal- 
ity condition. 

Solving for A instead of B has signifi- 
cant advantages, even though this involves 
taking another derivative. However, the 
total number of derivatives taken in the 
code is essentially the same. In fact, when 
centered finite differences are employed, 
Alfvcn waves are better resolved when A is 
used, because then the system of equations 
for one- dimensional Alfven waves in the 
presence of a uniform B^o field in a medium 
of constant density po reduces to 

Uz = (/ioPo)"^^xo^y, = B^oUz, (8) 

where a second derivative is taken only once 
(primes denote x-derivatives) . If, instead, 
one solves for the field, one has 

Uz = [ixoPoY^BmB'^, Bz = B^u'^, (9) 

where a first derivative is applied twice, 
which is far less accurate at small scales if 
a centered finite difference scheme is used. 
At the Nyquist frequency, for example, the 
first derivative is zero and applying an ad- 
ditional first derivative gives still zero. By 
contrast, taking a second derivative once 
gives of course not zero. The use of a stag- 
gered mesh would circumvent this difficulty. 
However, such an approach introduces ad- 
ditional complications which hamper the 
ease with which the code can be adapted to 
new problems. 

Another advantage of using A is that it 
is straightforward to evaluate the magnetic 



helicity, ( A • B) , which is a particularly im- 
portant quantity to monitor in connection 
with dynamo and reconnection problems. 

The main advantage of solving for A is 
of course that one does not need to worry 
about the solenoidality of the B-field, even 
though one may want to employ irregular 
meshes or complicated boundary condi- 
tions. 



4. Cache-efficient coding 

Unhke the CRAY computers that dom- 
inated supercomputing in the SOties and 
early 90ties, all modern computers have a 
cache that constitutes a significant bottle- 
neck for many codes. This is the case if large 
three-dimensional arrays are constantly 
used within each time step. The advantage 
of this way of coding is clearly the concep- 
tual simplicity of the code. A more cache- 
efficient way of coding is to calculate an 
entire timestep (or a corresponding substep 
in a three-stage 2N Runge-Kutta scheme) 
only along a one-dimensional pencil of data 
within the box. On Linux and Irix architec- 
tures, for example, this leads to a speed-up 
by 60%. An additional advantage is a dras- 
tic reduction in temporary storage that is 
needed for auxiliary variables within each 
time step. 



5. Large scale fields from helical turbulence 

Many astrophysical fiows are affected by 
rotation and gravitational stratification. 
These two effects can make the flows heli- 
cal. In the sun, for example, the kinetic he- 
licity of the flow is negative in the northern 
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Fig. 3. Visualization of the magnetic field in a 
three-dimensional simulation of helically forced turbulence. 
The turbulent magnetic field is modulated by a slowly 
varying component that is force free. 

hemisphere and positive in the southern. 
It has long been known that this can lead 
to the production of large scale magnetic 
fields §. 

Recent simulations of helically forced tur- 
bulence have shown that large scale fields 
are indeed produced 0] . These large scale 
fields are approximately force-free and of 
Beltrami type; see Fig. |. The prototype of 
a Beltrami field is B oc (cos 2;, sin 0), and 
it is easy to see that for this field J x B = 0, 
where J = V x B//io is the current density. 

At large scales, force-free Beltrami fields 
are only possible in a periodic box, as con- 
sidered in Ref. When the box in non- 
periodic, the large scale field can still be 
'nearly periodic'. In simulations with so- 
called 'pseudo-vacuum' boundary condi- 
tions, for example, a large scale field of the 
form 

B = (cos ^2;, sin 1^;, 0) cos 1^; (10) 

appeared; see Ref. ^ . We now consider the 
case of perfectly conducting boundaries. 
Unlike the case of pseudo- vacuum bound- 
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Fig. 4. Evolution of magnetic and kinetic energies, Em and 
Eki respectively, in a simulation with perfectly conducting 
boundaries. 

ary conditions, where the energy of the 
large scale field was somewhat below the 
kinetic energy, with perfectly conducting 
boundaries the energy of the large scale 
field can strongly exceed the kinetic energy 
of the turbulence; see Fig. 

The resulting super-equipartition of large 
scale magnetic energy is primarily a conse- 
quence of the fact that the large scale field 
is force-free and does not act back on the 
fiow. Force-free fields are generally helical 
and have to obey the equation of magnetic 
helicity conservation. 



d(A- B)/dt = -2r]fio{3 ■ B). 



'Ill 



The precise saturation level is obtained 
from the fact that in the steady state the 
magnetic helicity, (A • B), is constant. This 
implies that the current helicity, (J • B), 
vanishes. Now, splitting the field into large 
and small scale components, B = B + b, 
we have (J ■ B) = — (j ■ b). Assuming that 
the fields are fully helical at small and large 
scales, we have /io(J ■ B) = =i=A;i(B ) and 
/^o(j • b) = ±/cf(b^), where ki is the small- 
est wavenumber in the box and fcf is the 
wavenumber of the forcing. (Upper/lower 
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signs apply to positive /negative sign of the 
helicity of the forcing function.) Due to 
non-periodic boundary condition, however, 
the large scale field cannot be completely 
force-free. Therefore we allow for an effi- 
ciency factor eLs < 1 for the large scales 
(LS). Thus, the final balance equation is 



eLsfci(B') = fcf(b^), 



(12) 



i.e. the energy of the large scale field ex- 
ceeds that of the small scale field by a fac- 
tor eLgfcf/fci. This is consistent with Fig. ^. 
Note, for comparison, that for fully periodic 
boxes, eLs = 1 and therefore the factor of 
superequipartion is only kf/ki. 

The reason why the magnetic energy sat- 
urates in stages is connected with the oc- 
currence of different patterns of the large 
scale field: at early times the large scale field 
is dominated by a pattern with a relatively 
large wavenumber in one of the two hori- 
zontal directions (the boundaries are in the 
vertical direction). Initially the large scale 
fields has 8 nodal planes, but then it reduces 
to 4 and finally to 2 nodal planes. 

Finally, it should be emphasized that we 
have discussed here only the details of the 
saturation behavior. However, when the 
field is weak the magnetic field grows al- 
ways exponentially on a dynamical time 
scale (usually over many orders of mag- 
nitude), and is independent of magnetic 
helicity conservation. 



tral schemes. Explicit meshpoint schemes 
can readily be implemented on massively 
parallel architectures using High Perfor- 
mance Fortran (HPF) or the Message Pass- 
ing Interface (MPI). The 2A^-schemes of 
Williamson |I| are ideal for reducing the 
amount of storage while still allowing the 
temporal order of the scheme to be high. 
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6. Conclusions 



The use of high-order schemes proved to 
be a useful compromise between the cheap, 
but less accurate low-order methods and 
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